#ifndef Algebra_ODESolver_h
#define Algebra_ODESolver_h

#include "AlgebraModule.h"
#include "Object.h"
#include "ODESystem.h"

SG_NAMESPACE_OPEN

namespace Algebra
{
    // Forward declarations
    class ODESystem;
    class SecondOrderODESystem;

    /**
     * \addtogroup ODE
     * @{
     */
    /**
     * For structural dynamic problems, the following semi-discretized equations will be obtained
     * after the regular FEM operations in the space domain,
     * \f[
     *  \mathbf{M} \frac{d^2(u)}{dt^2} + \mathbf{C} \frac{d(u)}{dt} + \mathbf{K} u = \mathbf{f}
     * \f]
     *
     * What's more, the above ODEs can be transformed to a system of first order ODEs by
     * introducing the new variable \f$ v = \frac{du}{dt} \f$ and \f$ x = \left ( u,v \right )^{T} \f$.
     * And then the general motion equations can be re-written as followings
     * \f[
     *  \begin{matrix}I&0\\C&M\end{matrix} \begin{matrix} \frac{d(u)}{dt}\\ \frac{d(v)}{dt}\end{matrix} =
     *  \begin{matrix}f\\0\end{matrix} -\begin{matrix}K&0\\0&I\end{matrix} \begin{matrix} u \\ v \end{matrix}
     * \f]
     *
     * In summary, to solve time dependent problems assumes, the problem should be re-written in
     * the following form.
     * \f[
     *  \mathbf{F} \left(t, u, \frac{d(u)}{dt} \right)  = \mathbf{G} \left(t, u \right)
     * \f]
     *
     * After that, one ODE solver should be taken to solve for the solutions of the ODEs system.
     * During the long past time, a lot of third-party libraries in this sectors have been developed,
     * such as PETSc TS, Boost Odeint and so on.
     *
     * @see https://petsc.org/release/manual/ts/
     * @see https://www.boost.org/doc/libs/1_85_0/libs/numeric/odeint/doc/html/index.html
     *
     * @author nene
     * @date November, 2024.
     *
     */

    /** Basic Options for time integrators. */
    class ALGEBRA_EXPORT ODESolverBase : public Base::Object
    {
    public:
      explicit ODESolverBase ();
      virtual ~ODESolverBase ();

      /// Reset the current time.
      void SetTime (Real_t t);

      /// Set the maximum (or final) time for time stepping.
      void SetMaxTime (Real_t t);

      /// Reset the timestep at any time.
      void SetTimeStep (Real_t dt);
      /// Determine the current timestep.
      Real_t GetTimeStep ();

    protected:
      Real_t mCurrentTime;  //< the current time
      Real_t mFinalTime;  //< final time to step to
      Real_t mTimeStep;  //<  the current timestep size
    };

    /** Define the interfaces to solve a system of general ordinary differential equations. */
    class ALGEBRA_EXPORT ODESolver : public ODESolverBase
    {
    public:
      explicit ODESolver ();
      virtual ~ODESolver ();

      /// initialize the ODE solver.
      virtual void Initialize (ODESystem *odes);

      /// Perform a single time step from time t to time t+dt.
      virtual void Step (Vec_t &u);

    protected:
      ODESystem *mODEs = nullptr;  //< the representations of the ODEs
      bool       mIsInitialized = false;
    };

    /** Define the interfaces to solve a system of second-order ordinary differential equations. */
    class ALGEBRA_EXPORT SecondOrderODESolver : public ODESolverBase
    {
    public:
      explicit SecondOrderODESolver ();
      virtual ~SecondOrderODESolver ();

      /// initialize the ODE solver.
      virtual void Initialize (SecondOrderODESystem *odes);

      // /// time stepper for the specified ODEs M*d2udt + C * dudt + K * U = F
      // virtual void Step (const SpMat_t &M, const SpMat_t &C, const SpMat_t &K, const Vec_t &F, Vec_t &U, Vec_t &V, Vec_t &A);

      /// Perform a single time step from time t to time t+dt.
      virtual void Step (Vec_t &u, Vec_t &v, Vec_t &a);

    protected:
      SecondOrderODESystem *mODEs = nullptr;  //< the representations of the second-order ODEs
      bool                  mIsInitialized = false;
    };

    // /// Utility to register the ODE solver.
    // template <class CLASS>
    // class ODESolverProducer : public Base::Producer
    // {
    //   public:
    //     explicit ODESolverProducer (const char *className)
    //     {
    //         Base::Factory::GetInstance ()->AddProducer (className, this);
    //     }
    //     ~ODESolverProducer () override
    //     {}
    //     Base::Object *Produce () const override
    //     {
    //         return (new CLASS);
    //     }
    // };

    /** @} */
}  // namespace Algebra

SG_NAMESPACE_CLOSE

#endif  // Algebra_ODESolver_h